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Abstract: Long range rapidity correlations in A+A collisions are sensitive to strong color 
field dynamics at early times after the collision. These can be computed in a factorization 
formalism [1] which expresses the n-gluon inclusive spectrum at arbitrary rapidity separa- 
tions in terms of the multi-parton correlations in the nuclear wavefunctions. This formalism 
includes all radiative and rescattering contributions, to leading accuracy in a s AY, where 
AY is the rapidity separation between either one of the measured gluons and a projectile, or 
between the measured gluons themselves. In this paper, we use a mean field approximation 
for the evolution of the nuclear wavefunctions to obtain a compact result for inclusive two 
gluon correlations in terms of the unintegrated gluon distributions in the nuclear projec- 
tiles. The unintegrated gluon distributions satisfy the Balitsky-Kovchegov equation, which 
we solve with running coupling and with initial conditions constrained by existing data on 
electron-nucleus collisions. Our results are valid for arbitrary rapidity separations between 
measured gluons having transverse momenta p±, q± > Q s , where Q s is the saturation scale 
in the nuclear wavefunctions. We compare our results to data on long range rapidity cor- 
relations observed in the near-side ridge at RHIC and make predictions for similar long 
range rapidity correlations at the LHC. 
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1. Introduction 

In a high energy heavy ion collision, several thousand particles are produced in the initial 
interaction. The formation and evolution of the resulting fireball can be described in a 
framework where the incoming nuclei are sheets of strongly correlated coherent gluonic 
fields called Color Glass Condensates (CGC) [2-12], which are shattered in the collision 
to form strong classical fields called the Glasma [13-15]. The Glasma expands and ther- 
malizes to form a nearly perfect quark-gluon fluid, which eventually hadronizes and freezes 
out to produce the large observed multiplicity of particles. While there is a fair amount 
of circumstantial evidence on the temporal evolution of latter stages of this space-time 
scenario, at present it is the earliest times, with the strongest "Glasma" fields, that are 
most amenable to a systematic theoretical treatment. This is because the early time dy- 
namics at times of order 1/Q B < 1 fm is controlled by the saturation scale Q s , which is 
the characteristic momentum scale in the evolution of the bulk matter produced in the 
collisions [16, 17]. Estimates for the magnitude of are 1 — 1.4 GeV 2 for gold nuclei 
at RHIC and 2.6 — 4 GeV 2 for lead nuclei at LHC [18]. The existence of this semi-hard 
scale suggests that the Glasma may be described in weak coupling, thereby opening a new 
window into the study of strongly correlated quark-gluon matter. 

The properties of the Glasma can be investigated by measuring long range rapidity 
correlations of particles produced in the collision. This is because the requirement that 



- 1 - 



correlations be causal requires the latest proper time tj that two particles could have been 
correlated to be 1 



where the freezeout time Tf. Q . is the proper time at which particles from the fireball have 
no further interactions and Ay is the rapidity separation between the two particles. Thus, 
f° r r f.o. ~ 10 fm, two particles separated by 4 units in rapidity must have been correlated 
at no later than 1.4 fm. Strong correlations at space-time rapidity separations of Arf < 4 
units have been observed in the "near-side ridge" correlations measured by the PHOBOS 
experiment at RHIC [19]. Correlations up to A77 = 1.5 have been extensively studied by 
the STAR collaboration [20]. At the LHC, multi-particle correlations at very large rapidity 
separations can be studied; these are correlations that must have been created at proper 
times well below a fermi, thereby providing a unique window into the non-linear dynamics 
of strong classical color fields in QCD. 

It is therefore interesting to study the nature of these correlations and what they reveal 
about the Glasma. Further, since these correlations occur at very early times, after the 
collision, they are closely related to multi-parton correlations in the nuclear wavefunctions 
themselves. Multi-particle production in the Glasma, and its relation to correlations in 
the nuclear wavefunctions can be described in a weak coupling QCD framework where the 
degrees of freedom are strong color sources p a ~ 1/g in the nuclei (where g is the QCD 
coupling constant) and gauge fields. Before the collision, the distribution of sources and 
fields in the nuclear wavefunctions evolves with rapidity; the evolution equations for the 
color source distributions are the JIMWLK renormalization group equations [5-12]. After 
the collision, the color sources become time dependent, thereby enabling particle production 
in their radiation field. In Refs. [21,22], a field theory formalism was developed to compute 
moments of the multiplicity distribution in the Glasma systematically as an expansion 
in powers of g 2 , while simultaneously resumming contributions of order gp ~ 0(1) from 
arbitrary numbers of insertions of color sources at each order in g 2 . 

The naive expansion in powers of g 2 however breaks down because at each order there 
are large logarithmic contributions in x\ t 2, the momentum fractions of partons in each of 
the nuclei, such that g 2 h^l/xi^) ~ 0(1) at small x\^- These contributions therefore have 
to be resummed as well. In Refs. [1,23,24], it was shown that inclusive observables 2 in the 
Glasma can be expressed in a factorized form 



are strongly correlated. 

2 This factorization is proven, to leading logarithmic accuracy in x, for inclusive multi-gluon spectra. It 
is straightforward to check that it applies to the expectation value of local operators such as the energy- 
momentum tensor T M ", and multi-point correlations of such operators. It is unlikely to apply to exclusive 
final states that impose a veto on particle production in some regions of phase-space. 




(1.1) 




lr rhis expression is valid in the scenario when the space-time rapidity and momentum space rapidities 
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where the Wilson lines 

/•''.; _ r i 

^1,2(2/, x_l) = P expig / dz T — T p 1)2 (z T ,x±) . (1.3) 
Jo v ± 

are ordered in rapidity (or, equivalently, in the longitudinal coordinates x T ). In this defi- 
nition of the Wilson line, the rapidity y is measured from the fragmentation region of the 
projectiles 3 . The p% t 2 are the color source densities of the nuclei in Lorenz gauge at a given 
transverse co-ordinate xj_ and longitudinal position . The W J s are universal weight 
functionals (diagonal elements of density matrices) that give the probability distribution 
of a given configuration of sources (or equivalently the Wilson lines O]^)- 

If the separation scale between the sources and the fields in one of the nuclei (moving 
in the + direction) is A + , the requirement that the physics be independent of this cutoff 
gives rise to the JIMWLK renormalization group equation [5-12] 

w A+ [n 1 ] = -n A+ w A+ [n 1 ], (1.4) 



a In (A 

where is the JIMWLK Hamiltonian at the cutoff scale A + . The precise form of the 
JIMWLK Hamiltonian is not needed in the rest of this paper, and we refer the reader to 
Refs. [5-12,23] for further details. The distributions Wf^i^] that enter in eq. ( |1.2j ) are the 
limits when A± — > of the cutoff dependent distributions W\ ± [fij.,2] 



The "master" formula in eq. (1.2) is valid for all moments of the multiplicity distribu- 
tion of gluons produced in the Glasma. Given a non-perturbative initial condition, the W 
weight functionals encode information about gluon correlations at all transverse positions 
and rapidities. A remarkable feature of eq. ( |1.2p is that the only "process dependent" input 
on the right hand side of the expression is the observable computed to leading order 4 in a s ; 
O hQ , albeit non-perturbative, is obtained by solving classical Yang-Mills equations for the 
two nuclei and has been studied extensively for the single inclusive [25-33] gluon spectra 
using numerical lattice methods. 

The Glasma Flux Tube picture of A+A collisions [34] is a consequence of this master 
formula. At short times after the collision, the solutions of the Yang-Mills equations 



3 It is therefore different from the usual laboratory frame rapidity y used as a measure of the longitudinal 
momentum of a particle in the final state of a reaction. For projectile 1, moving in the +z direction, they 
are related by y — Ibeam — y, and for projectile 2 by y = — Ybeam — V- The rapidity difference from the beam 
y is related to the upper bound in eq. ( |l.3| ) by y — ln(P ± x :F ) where P ± denotes the total longitudinal 
momentum of the projectiles 1 and 2 respectively. 

4 By leading order, we mean the first term in the expansion 

° H = [ Co + Cl # 2 + C2 9 4 H ] . ( L5 ) 

where each term corresponds to a different loop order. Each of the coefficients c n is itself an infinite series 
of terms involving arbitrary orders in {gpi,2) v - The "leading order" contribution, 

°lo[P1,H = JL: (1-6) 

corresponds to an infinite sum of tree diagrams. 
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provide longitudinal chromo-electric and chromo-magnetic fields in the forward light cone. 
In the McLerran-Venugopalan (MV) model [2-4], the distribution of color correlations in 
the nuclear weight functionals W is Gaussian with its width set by the saturation scale Q s . 
The averaging over the weight functionals in eq. ( |1.2| ) ensures that only chromo-electric 
and magnetic fields localized in transverse areas of size 1/Qg contribute to multi-particle 
production obtained by averaging over all events. In the MV model, color correlations 
extend to arbitrarily large rapidity intervals; this results in a picture of multi-particle 
production as arising from boost invariant flux tubes of size This picture, albeit 

simple, gives a qualitative [34] and even semi-quantitative [35,36] description of the near- 
side ridge correlations observed in the RHIC heavy ion experiments. 

Reality however is more complex and the boost invariance of the Glasma flux tubes is 
violated both by quantum evolution effects (real gluon emissions and virtual corrections) in 
rapidity between the beam rapidity and the rapidities of the measured gluons and likewise 
by quantum evolution between the measured gluons. When the maximum rapidity interval 
between measured gluons AY <C l/a s , quantum radiation between the observed gluons is 
not significant and correlated gluon emission is approximately independent of AY. The 
factorization formalism for this case was developed in Refs. [23,24]. A quantitative un- 
derstanding of how correlations depend on Ay, for arbitrary Ay < 2Yb eam , requires that 
one understands the dynamics of real and virtual quantum corrections between the tagged 
gluons. As demonstrated in Ref. [1], all the necessary information is contained in eq. ( |1.2j ) 
and the general formula for two gluon correlations was derived in that paper. 

In this paper, we will exploit the formalism of Ref. [1] to evaluate the rapidity de- 
pendence of correlated two gluon emission in A+A collisions at RHIC and LHC energies. 
Because solving the JIMWLK equation is highly computationally intensive [37], we will 
instead use a mean field approximation to this evolution equation known as the Balitsky- 
Kovchegov (BK) equation [38,39]. The BK equation is a very good approximation to the 
JIMWLK equation [40], albeit it must be noted that this may be true only for a limited 
class of observables. Recently, significant progress has been made in computing Next-to- 
Leading-Order (NLO) contributions to the BK equation [41,42] and the results have been 
successfully applied to phenomenological studies of HERA DIS data at small x [43,44]. 
To apply this framework to nuclear collisions, we will first fix the initial conditions for the 
running coupling BK evolution of unintegrated gluon distributions for nuclei by fitting the 
existing inclusive e+A fixed target data. We will then apply our results to compute the 
rapidity dependence of two gluon correlations in A+A collisions. 

Our paper is organized as follows. Readers interested primarily in the results of the 
paper can proceed directly to section ||. In section |2|, we will restate the results of Ref. [1] 
for the single and double inclusive gluon spectrum and demonstrate how the expressions 
simplify vastly in the mean field BK approximation of quantum evolution. In section |3], the 
correlated two gluon distribution is expressed in terms of the unintegrated gluon distribu- 
tions of nuclei. These distributions are determined in section ||| from numerical solutions of 
the BK equation with running coupling. To constrain the initial conditions for the evolu- 
tion of these nuclear unintegrated distributions, we fit available fixed target e+A data - 
the initial conditions for proton unintegrated distributions were determined previously in 
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global fits to the HERA e+p data [43,44]. In section [|, we discuss our results for correlated 
two gluon production in the context of RHIC data from STAR and PHOBOS. We make 
predictions for what one may expect at the LHC. The final section contains our conclu- 
sions. Details of the computations, solutions of running coupling BK equations and fits to 
e+A fixed target data are contained in two appendices. 



2. Double inclusive gluon spectra at arbitrary rapidities 

The general leading log factorization formula (|h^) gives the distribution of gluons at all 
rapidities in the problem [1,45,46], including all the rapidity correlations in the leading 
log approximation. For single and double inclusive gluon production we only need the 
distribution of Wilson lines at one or two gluon rapidities respectively. We shall now 
specialize the generic formula (1.2) to the this specific case. The single inclusive gluon 
spectrum dAq/d 3 p at leading order depends only on the Wilson lines 5 f2i i2 (y, x ±) at the 
rapidity y = y p of the produced gluon and not on the whole rapidity range contained in 
eq. ( |1.2[) . Therefore, we can simplify eq. (|1.2| ) by inserting the identity 



1 = J [DU 1>2 (x±)\ 5[U 1;2 (xi_) - fii, 2 (lfc,xi)] (2.1) 

and by defining the corresponding probability distributions for configurations of Wilson 
lines at the rapidity y p 

Zy.px^)} = J [£>n lj2 (y,x ± )] W[n lt2 (y,x x )] 5[Ui, 2 (x ± ) - 1)2 (y p ,x ± )] . (2.2) 

One then obtains the all order leading log result for the single inclusive gluon spectrum at 
the rapidity y p to be 



dNl ^ - f [DU lM DU 2 (x ± )) Z yp [UJ Z yp [U 2 ] dNl ^ 



(2.3) 



^p ± dy p /^ og J l — y — -v—yj i-J ~» i-J d 2 p±dyp 

Note that the distribution Z y [U] obeys the JIMWLK equation, 

dy p Z yp [U]=H yp Z yp [U], (2.4) 

which must be supplemented by an initial condition at a rapidity close to the fragmentation 
region of the projectiles. Eq. 2^3 is illustrated in figure |]. 

As will become clear shortly, it is more convenient to express eq. ( [2 .3D in terms of color 
charge densities as 



< = / m w*- ji z " w z " w 



(2.5) 

LO 



J From here onwards, we shall use the standard definition of the rapidity y instead of the rapidity distance 
y from the projectile. In terms of the Wilson lines introduced in the previous section, this translates to 
Qi(2/,x_l) = Qi(Ybeam - V,x-±-) and Q 2 (j/,xj_) = Q.i{y + lbcam,x±). 
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Initial configuration 
for nucleus 1 



JIMWLK evolution 
for nucleus 1 



[dNx/d 3 ?], 



Figure 1: Diagrammatic representation of the various building blocks in the factorized formula 
for the inclusive single gluon spectrum. The lower part of the figure, representing nucleus 2, is made 
up of identical building blocks. 



Here, the superscript on p\ 2 denotes the color charge distribution of nucleus 1 or 2 evaluated 
at the rapidity y p . 

Following a similar though slightly more involved derivation [1], one obtains from 
eq. ( |1.2D the expression 

dN 2 



d 2 p±dy p d 2 q±dy q 



[Dpl(x ± )Dp p 2 (x ± )Dpl(x ± )Dpl(x ± )} 

x z v v G v P ,y q K> Pil z v q IpD G y q ,y P \p%(%\ 
dJVi[p?,^]| 



X 



d 2 p±dy p 



d 2 q ± dy q 



(2.6) 



It is important to note that here we have taken y q > y p , where y p is at an earlier stage 
in the evolution of projectile 1 and likewise, y q is at an earlier stage in the evolution of 
projectile 2. This convention will be followed for the rest of this paper. Also, in eq. (p.q*), 
G y q ,y P [pi 2' Pi 2] i s a Green's function of the operator d y — Tt y , 

dy q G yq ,y v [U q , U^] = H yq G yq , yp [17* U?] , (2.7) 

with the boundary condition 

I™ G Vq>Vp [p^(P)=8[(fl-(P] . (2.8) 
y q ~ *yp 

This Green's function describes quantum evolution between two specified rapidities, in 
the presence of strong color fields from the projectiles. It relates the distribution of color 
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sources at a given rapidity with the distribution at another rapidity through the relation 



Z yp \& = J [Dp q 2 ] Z yq [pl\ G yq!Vp 



(2.9) 



At this stage, it is important to note that the BK mean field form of the Balitsky- 
JIMWLK evolution equation for the two point Wilson line correlator (the "dipole cross- 
section" ) requires that the correlator of the product of traces of two pairs of Wilson lines 
factorizes into the product of the correlators of traces of pairs of Wilson lines 6 , 

tr[E^(xj_)£f(z_L)] tr[C/t(z ± My ± )] 

tr[f/t( x± )L>(z ± )]) (tr[^(z ± )Lf(yj.)]) , (2.10) 

to leading order in a 1/N C expansion. The averages ( • • • ) are performed over the color 
sources of a large nucleus with the weight functional Z y . The factorization in eq. ( 2.1Q ) 
can be achieved with Gaussian correlations among the color sources. Note however that we 
need a non-local Gaussian distribution to accommodate the quantum BK evolution [47]. 
One obtains therefore 7 



Z, 



Up 



exp 



P°2( x -L)Pl,2(y-L) 

,y± /4 12 (%>,x ± -y±) 



(2.11) 



where p\ 12 (y p ,x± — y±) represents the color charge squared per unit area of nucleus 1 
or nucleus 2 as seen by a particle having rapidity y p . Even though in this work we will 
consider collisions between identical nuclei, we shall retain the explicit notation for 
generality. 

From eq. (2.9), because the Z y functionals on the l.h.s and the r.h.s are both Gaussians, 



the Green's function G yqt y p must be Gaussian as well. One obtains 



°y P ,y q bi.Pi] = ex P 



G y g ,y P \pb(%\ = ex P 



A Pl (x ± )A Pl (y ± ) 
2 7x ±>yx A/i^(xj. -y ± ) 

Ap 2 (x ± )Ap 2 (y ± ) 
AmL(x± -y±) 



(2.12) 



where we have defined 



Api(x ± ) = pl(x±) - P?(x_l) 

A/9 2 (XJ_) = /^(X ± ) - ^(X ± ) 

(t± ) = fi 2 Al (Vg ) r -L ) - Pm (Vp » r ± ) 
(rj.) = p%(y P ,r ± ) - /4 2 (y<z,r±) 



(2.13) 



6 These Wilson lines are defined in terms of the color charge densities through eq. (1.3) projected on to 
a particular rapidity U(x±) = exp lig-^ r p a (x±)t a ) , where the t a, s are the generators of the fundamental 



representation of SU(N C ). The corresponding expression in the adjoint representation is given in eq. (3.6) 
7 It is to be understood that repeated color indices a are summed over. 
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Note that because of our choice y q > y p , Afi 2 as defined is always positive 8 



Because the Green's functions in eq. (2.12) are expressed naturally as Gaussians in the 



new variables introduced in eq. ( |2.13| ) , we can rewrite our general expression for the double 
inclusive distribution as 



dN 2 



d 2 p±dy p d 2 q_Ldy g 



[Dp?(x ± )£p«(x ± )£Ap 1 (x ± )DAp 2 (x. L )] 



x z y P iPi] G v v ,y q [Pi> Pi] Z v q [Pl\ G y q ,y P \p\ > f%\ 



d 2 p±dy p 



WW,/® 



LO 



d 2 q±dy q 



(2.14) 



With the ZyS from eq. ( 2.11j ) and the G yq>y 's from eq. ( |2.12| ), the only ingredient missing 
in obtaining a final expression for eq. ( |2.5| ) and eq. ( |2.14| ) is the expression of the leading 
order single inclusive spectrum in terms of the color charge densities of the two nuclei. This 
expression and the subsequent simplification of our equations for the inclusive distributions 
will be discussed in the next section. 

Initial configuration 
for nucleus 1 

JIMWLK evolution 
— > for nucleus 1 
from Y beam to Y p 



[dNxAPp] 




JIMWLK evolution 
for nucleus 1 
from Y n to Y n 



[dNx/d'q], 



Figure 2: Diagrammatic representation of the various building blocks in the factorized formula 
for the inclusive 2-gluon spectrum. As in the previous figure, the corresponding evolution from 
nucleus 2 at the bottom of the figure is not shown explicitly. 



3. Gluon correlations from unintegrated gluon distributions 

The leading order single particle inclusive distribution, for a fixed distribution of sources, 

8 /j,\ is proportional to the saturation scale, and therefore increases as one evolves away from the frag- 
mentation region of a projectile. 
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is given by 

dNi [p!,p 2 ] 



d 2 p±dy p 



1 lim f d 3 xd 3 y e^*"") (Sg - i£ p )(a° + i^) 



x ^(p^p) A*(x)[ Pl ,p 2 ) A*(y)[ Pl ,p 2 ] . (3.1) 

A. a 



p*AZ( P ) = - lfabc g* I ^iAp,k ± r ^ m ^ (3.2) 



The gauge fields A®(x)\pi, p 2 ] are solutions of the classical Yang- Mills equations in the 
forward light cone after the nuclear collision for a fixed configuration of sources pf 2 i n each 
of the nuclei. For Fourier modes kj_ of the color charge densities which obey pi 2(kj_)/k^ > 
1 (which is the case for Q s > k±), only numerical solutions for A^(x) are known [25-33]. 
However, for pi^/k^ <C 1, valid for Q s <^ k±, one can perturbatively expand the gauge 
field in powers of pi t2 /k 2 ± and one obtains [48,49] 

(2tt)2 VP ' ^ ki( P± -kj 

where f a \, c are the SU(3) structure constants and pi 2 are the Fourier transforms of the 
color charge densities of the two nuclei. Here is the well-known [50,51] Lipatov vertex 9 . 

For the single inclusive distribution in eq. (|2.5| ), using eq. ( |3. 1[ ) and eq. ( |3.2| ) and the 
correlator 

(p a (k x )f?(k' ± )) = (2irfp\(y) 5 ab 5(k ± - k' ± ) , (3.3) 

one obtains 

diVi \ c 2g 6 iV e (jV c 2 - 1) 1 f p 2 Al (y P ,k ± )p 2 A2 (y P ,P±-k ± ) 

~ b ± TnZvs ZT / d k ^ T7TTZ H v? > ( 3 - 4 ) 



d2p ± dy p / LLog " (2tt)5 pi7 " ki( Pl 

where S± is the transverse area of the overlap between the two nuclei. The unintegrated 
gluon distribution can be expressed as [47,52,53] 

K^ k ^ = / d ' x± £ik " x± ( Tr ( c/t (°) c/ ^))) ' < 3 - 5 ) 

where the matrices U are adjoint Wilson lines evaluated in the classical color field created 
by a given partonic configuration of the nuclei A\ or A 2 . For a nucleus moving in the — z 
direction, 



U (xj_ ) = P + exp 



+00 

ig f dz + ^ TPa (z + , x± )T a 



(3.6) 



Here the T a are the generators of the adjoint representation of SU (N c ) and P+ denotes 
path ordering along the z + axis. At large k^, the Wilson lines can be expanded in powers 
of the sources to give, for Gaussian correlations, 

4> A (3,M = M*#M - 1)%^ . (3.7) 



9 The components of the Lipatov four vector are Z/ + (p,kj_) = — L (p, k^) = ^-^ — — ^ — — 
i l (p,k ± ) = -2k l ± . 
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Substituting this relation in eq. fl3.4|) , we obtain the well known kj_ -factorization expres- 
sion [54,55] for the single inclusive gluon distribution valid for p± ^> Q s : 



dN x \ 2a s N c S± 1 f d z k 



J J2n^ ®A 1 (y P ,'k ± )<f>A 2 (yp,P± - k±) , (3.8) 



\d 2 P± dy p / LLog 2n*(N*-l)p*_J (2tt) 

where we denote &a = ^/(vri?^) to be the unintegrated gluon distribution per unit of 
transverse area. 

The corresponding expression for the double inclusive distribution is more involved. 
The r.h.s of eq. ( |2.14 ) has the product of two single inclusive distributions, one for a 



gluon with three momentum p and likewise another for a gluon with three momentum 
q. From eq. ( |3.1| ), this corresponds to the product of four gauge fields. As for the single 
inclusive case, the double inclusive gluon spectrum can be computed numerically using 
lattice techniques where Yang-Mills equations are solved to obtain the gauge fields as a 
function of proper time after the collision. This computation has been carried out recently 
for the Gaussian MV model [56] . Because this model does not include the effects of small x 
evolution, it is not ideal for the purpose of investigating the dynamics of long range rapidity 
correlations. Incorporating small x evolution effects in the non-perturbative computation 
is outside the scope of the present work. We will instead consider here, as in the previous 
discussion of the single inclusive distribution-see eq. ( |3.2[ ) , the perturbative limit of p± , 
q± ^> Qs, where the gauge fields can be expanded as bilinear products of the color sources 
of the two nuclei. The dependence of the leading order double inclusive gluon spectrum on 
four gauge fields, then translates, in this perturbative limit to the product of eight color 
charge densities. The averages over color sources in eq. ( p. 14 ) are therefore averages over 
the general matrix element 

^e/^ (pjq;{k . ±}) s ^/.P (k2±) ^(^ ± )^(k 1± ) (k 3± ) 

X P*2 9,P (.P± ~ k 2-0 P2 ,q (l± ~ k 4_L)$>' P (P± - ki_L) P2' q {(l ± - k 3J J 

(3.9) 

where we denote by a superscript p or q the rapidity which the color sources correspond 
to. Further, these products of gauge fields contain bi-linear scalar products of the Lipatov 
vertices. These can be simplified [34, 57] and expressed as 



g(p,q;{k i± }) 



16 [(ku • p ± - k 2 1± ) (k 2 _L • p± - \4 ± ) + (leu. x p ± ) • (k 2 ± x p±)} 



(2^) 8 k2 ±k 2 ±P i(p ± -k 1± ) 2 ( P± -k 2± ) 2 

[(k 3 _L • q_L - k|_J (k 4± • q ± - k| ± ) + (k 3± x q ± ) • (k 4± x q ± )] 

kL^qK^-L - k 3±) 2 (q± - k 4 ±) 2 

(3.10) 
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The double inclusive distribution in eq. ( 2.14| ), for transverse momenta p±-,q± 3> Q s , 
can therefore be expressed as 



dN 2 \ g 



,12 



jabc ja'de j a $9 fofhi 



d 2 P±dy p d 2 q ± d y(? / TT _ 16(2tt)6 

r 4 

j Q(p,^{^±})^ bcdef9M (p,^,{^±}) , (3.H) 



x 



i=l 



in terms of jF 6cde ^ hl (p, q; {kjj_}) and G(p, q; {kjj_}) defined above. 

We shall now sketch how one evaluates these quantities with further details of the 
computation given in appendix A. We begin with the evaluation of the color averages in 
jrbcdefghi^ q. {kjj.}). Because the Z's in eq. ( |2.11| ) and the G in eq. fl2.12| ) are Gaussian 
weight functionals, the relevant color source correlators are the equal rapidity correlators 



^(kjjp^)) = (2^) 2 rt 2 (k ± - k'J^^kj 
^ a ' 9 (k ± )p*/(k' ± )\ = (2vr)W(k ± - k! L )n\{y q ,k A 



(A^ 2 (k ± )A^ i2 (k' ± )> = (2vr) 2 rt 2 (k ± -k' ± )A M i 12 (k ± ) , (3.12) 



,10 



and the non-equal rapidity correlators 

;^(k ± )p^( k ' ± )) = ((ApUk ± ) +pr p (k ± )) 

= (2^ 2 ^5 2 (k ± -k' ± ) / u 2 4i (y p ,k ± ) 

;^(k ± )^(k' ± )) = (n a >«(k ± ) (Ap 2 (kD + p 6 /(k' J 

The correlators for the dependent variables p\ and (F 2 ( see e Q- (|2.13|) ) are 

>r ,<? (k±)/^(kl)) = (2vr) 2 rt 2 (k ± - k' ± ) ^ Al (y q ,k ± ) , 
y 2 a ' p (k ± )f% p (k' ± )) = (2vr) 2 rt 2 (k ± - k' ± )» A Jy p ,k ± ) . (3.14) 



(2vr) 2 rt 2 (k ± -k' ± )^ 2 A (y q ,k ± ) . (3.13) 



With the relations listed, we can now evaluate eq. ( |3.9|) . Simple combinatorics gives 
us a total of 9 possible pairwise contractions in eq. ( [3.9D , The details of these are listed 
in appendix A. One of the contributions (eq. (A.l)) gives the non-correlated contribution 
to the two gluon inclusive distribution. Subtracting this term therefore results in the 
correlated two gluon inclusive distribution 

c(p ' q) = (pS^i} - (d^pr) (d^qi) • (3 - 15) 

When one evaluates the other 8 terms that contribute to C(p, q), one observes that only 4 of 
these give leading contributions. The <5-function contributions from these terms (eqs. ( |A.2j ), 



3 This follows from the vanishing of terms odd in p or Ap. 
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( |A,3|) , ( |A.4j) , (| A . 5|) ) give kij_ = k 2 j_ and k 3 j_ = k 4 j_. Substituting this into the expression 
for G{p, q; {kjj_}) in eq. ( 3.1C| ), one finds that it simplifies considerably to read 



S(p,q;{k 4± }) 



16 



(27T)* kl ± kl ± p\ q l (p ± - k 1± ) 2 (q ± - k : 



(3.16) 



3±y 



Combining the four leading contributions from _7 rbc(fe /9' w (p ) q; {kjj_}) and the corresponding 
expressions from G(p, q; {kjj_}) to C(p, q) (see eqs. ( |A.10| ), ( |A.11[ ), ( |A.12[ ) and ( |A.13j ) in 
appendix A), we obtain 



C(p,q) 



aj Ni(NZ-l)S ± 



d 2 k 1± x 



16^° d\ piqi 

j$A 1 (yp ) ki_ L )$A 2 (y P ,P_L - ki ± ) [$A 2 (%,q± +kij.) + $A 2 (y g ,q± -kij.)] 

+®m (.y q , ki±)$Ai {y P , p± - ki±) [*ai q± + ki±) + $Ai (v g , <u - k i± )] | , 

(3.17) 

where the <3?'s are unintegrated gluon distributions per unit of transverse area and cIa = 
— 1. We have used here the relation between /x 2 and the unintegrated gluon distribution 
<3? given in eq. ( |3.7| ). This expression is the central result of this paper 11 . We have obtained 
an expression for the double inclusive gluon distribution, valid to all orders in perturbation 
theory to leading logarithmic accuracy in x and for momenta p±, q± 3> Q s , entirely in terms 
of the unintegrated gluon distributions of the two nuclei evaluated at the rapidities y p and y q 
where y p < y q . The corresponding expression for y p > y q is obtained by replacing A\ <-> A2 
and y Pi q — > —y p , q - We should emphasize that the notation used in eq. ( 3.17 ) stipulates that 
the un- integrated gluon distributions are evaluated at rapidities y PjQ ± Ybeam- 



4. Running coupling BK evolution 

In the previous section, we established that the correlated two gluon spectrum for arbitrary 
rapidities can be computed in terms of the unintegrated gluon distributions of the two 
nuclei evaluated at these rapidities. In this section, we shall discuss how one computes this 
unintegrated gluon distribution and its evolution with x. In the next section, we shall use 
the results for the unintegrated gluon distribution to evaluate eq. ( 3.17] ) for the correlated 
inclusive two gluon distribution. 



In eq. ( |3.5| ), we defined the unintegrated gluon distribution in a nucleus in terms of 
the correlator of two adjoint Wilson lines averaged over the color charge distribution in a 
nucleus. Because these averages ( • • • ) and those of correlators of fundamental Wilson lines 
are Gaussian correlators in the large limit, one can express these correlators respectively 



11 We note that expressions for the double inclusive cross-section have been previously derived [58] within 
the framework of Local Reggeon Field Theory [59]. At present the connection between the two frameworks 
is completely unclear. 
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as [52,60] 

Tr (U(0)tf(r ± )} = N 2 C e ~c A nr^Y) 



Tr (u(0)U\r ± )) = N c e -C F ?(r±X) ? (4.!) 



where C A = N c is the Casimir in the adjoint representation and C F = (iV 2 — 1)/2N C is 
the Casimir in the fundamental representation. The function T is closely related [47] to 
the variance of the non-local Gaussian weight functional in eq. ( |2.11 ) and is therefore the 



same in both the fundamental and adjoint cases. One can therefore, in the large N c limit, 
simply express the correlator of two adjoint Wilson lines as the square of the correlator of 
two fundamental Wilson correlators. 

The correlator of two Wilson lines in the fundamental representation is simply re- 
lated to the dipole amplitude for the scattering of a quark-antiquark dipole (of transverse 
separation rj_) off a nucleus as 12 

T(r ± ,Y) = 1 - i-Tr (u\0)U(r ± )^ , (4.2) 

where U is a Wilson line in the fundamental representation. (See our previous discussion 
of these in the context of eq. ( 2.10[) .) Using eq. (|4.l| ), one can write the unintegrated gluon 



distribution in the adjoint representation (per unit of transverse area) in eq. (|3.5| ) as 

*Ai, a (a?,fcj.) = / r ± dr ± J (k ± r ± ) [l - T Al 2 (r ± , ln(l/x))] 2 . (4.3) 



We therefore need to determine the dipole amplitude T and its evolution with rapidity 
Y (=\ii(xq/x)) as an input in eq. (|4.3j ) to extract the unintegrated gluon distribution. The 
dipole amplitude is obtained from the Balitsky-Kovchegov (BK) equation [38,39], which 
is a non-linear evolution equation describing both gluon emission and multiple scattering 
effects in the interaction of the quark-antiquark dipole with a nucleus in the large N c limit. 
It can be expressed as 



dT(r,Y) f 

-= / dri /C LO (r,ri,r 2 ) x 

[T(n, Y) + T(r 2 ,Y) - T(r, Y) - T(r u Y) T(r 2 , Y)] , (4.4) 



dY 



with the leading order BFKL kernel [61] given by 

/C LO (r,n,r 2 ) = ^-^| , (4.5) 

where r 2 = r — n. As we discussed previously, the BK equation for the amplitude is 
equivalent to the corresponding JIMWLK equation [5-12] of the Color Glass Condensate, in 
a mean field (large N c ) approximation where higher order dipole correlators are neglected. 



2 We assume translation invariance in the transverse plane to set the quark transverse coordinate to zero. 
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In the context of the BK equation, the leading order kernel corresponds to resumming 
the leading (a s ln(xo/x)) n terms arising at small x from all orders in perturbation theory. 
It is well known however that running coupling contributions qualitatively modify the small 
x evolution beyond leading logarithms in x and there has been considerable recent work 
to include these corrections to the BK equation [41,42]. The running coupling equation 
describing the evolution of the dipole amplitude however takes exactly the same form as 
eq. Q4.4D with a modified evolution kernel given by 

a,(r)iV e [ r 2 1 ( a,(4) ,\ , 1 f <>,(4) ,X\ B , 

The subscript in JC-q&i. refers to the "Balitsky prescription" for the evolution kernel, which 
corresponds to a scheme where some particular ultra-violet finite terms are also included 
along with the running coupling contributions to make the remainder numerically less 
important. For a more detailed discussion, we refer the reader to [62]. In this work, the 
NLO contributions not encompassed by the kernel in eq. ( |4.6| ) will be ignored. As argued 
previously [43], these contributions are systematically smaller than the running coupling 
contribution included here, especially at large rapidities. 




Figure 3: Unintegrated gluon distribution in the adjoint representation at Y = 0, 2, 6, 10, 15 (from 



left curve rightwards) with the Balitsky prescription for the kernel in eq. (4.6) as well as for the 
fixed coupling case. The distribution is in units of N c irR\/ 'a 8 . 

In fig. H|, we show results for the unintegrated gluon distribution versus transverse 
momentum squared determined from the evolution with rapidity of the dipole amplitude 
in the adjoint representation (see eq. fl4,3|)) with i) the fixed coupling BK kernel, and ii) 



with the Balitsky prescription for the kernel in eq. (4.6)). As we will describe below, the 
initial conditions for the latter figure are constrained by fixed target e+A data. We note 
that the evolution of the unintegrated gluon distribution with Balitsky's prescription for 
the running coupling effects is significantly slower than the evolution with a fixed coupling 
constant. 
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Figure 4: The x and Q 2 dependence of the normalized ratio of structure functions F^ in nuclei. 
The curves in the left figure includes effects due to the small x evolution of the dipole cross-section 
described by the BK evolution with the modified kernel in eq. (|t.6|) . The curve in the right figure 
is sensitive to the Q 2 dependence of the initial condition alone because it is evaluated at relatively 
large x. Details regarding the parameters of the initial condition are discussed in appendix B. The 
data are from the NMC collaboration [63] . 




Figure 5: The A dependence of the ratio of structure functions given by data from the NMC 
collaboration [64]. The corresponding curves for other initial conditions are in appendix B. 

The BK equation with the modified kernel in eq. fl4.6|) was first applied in Refs. [43,44] 
to a phenomenological study of the HERA data on the proton structure function F2. Two 
sets of initial conditions for the dipole amplitude at the initial rapidity Y = Yq were used- 
the GBW [65] and MV initial conditions [2-4]— and their parameters determined from fits 
to the HERA data. To constrain the initial conditions for nuclei and therefore extract 
the nuclear unintegrated gluon distribution, we performed a fit to the available NMC data 
on the nuclear structure function F^fe, Q 2 ). The details of the fit and the results are 
described in detail in appendix B. We show here in figs. |4] and |5] representative plots of fits 
to x, Q 2 and A dependence of the fixed target e+A data. Good fits to the available data 
are obtained for both sets of initial conditions for particular parameters. With the initial 
conditions for the BK equation fixed by the NMC data, we shall now use the corresponding 
unintegrated gluon distribution to study long range rapidity correlations in the Glasma. 
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5. Results for long range rapidity correlations in the Glasma 



In this section, we will make use our result in eq. ( 3.17] ) for the double inclusive gluon 



distribution to compute long range rapidity correlations in A+A collisions at RHIC and 
the LHC. The essential ingredient in eq. ( p. 17 ) is the unintegrated gluon distribution which, 



as shown in eq. (4.3), is simply related to the dipole scattering amplitude. The evolution 



of the dipole scattering amplitude with rapidity (or x) is described by the BK evolution 



equation given in eq. (4.4), with the modified kernel given in eq. (4.6). The rapidity 
dependence of the double inclusive gluon spectrum therefore provides a sensitive test of 
high energy QCD evolution. 



Equation ( |3.17 ) is derived in the leading \nx approximation, where all transverse 



momenta are assumed to be parametrically of the same order as Q s . In this approximation 
the x-values at which the unintegrated gluon distributions are evaluated are not exactly 
determined, as long as x ~ e ±y Q s / ' \fs, where y is the appropriate rapidity of the produced 
gluon (y p or y q ) and the sign depends on the nucleus (1 or 2) considered. We define the 
longitudinal momentum fractions of the produced gluons with respect to nucleus 1 or 2 
(denoted by subscripts) 

T1 — E^p-yp ■ Tl — !LL P -y<i 

x 2p = *±e** ; x 2q = ^e + y (5.1) 

In the above expression, p± and q± are the transverse momenta of the produced gluons. 
The unintegrated gluon distributions with momentum argument pj_ ± ku_ and qj_ ± ku_ 
in eq. ( 3.17j ) are evaluated at these values of the momentum fraction. For the uninte- 



grated distribution with momentum argument ku_ we replace the transverse momentum 
in eq. ( |5.1| ) by (pj_ + q_i_)/2 to make our evaluation of eq. ( |3.17 ) manifestly symmetric in 



p_i_ and q_i_ 13 . Our derivation in Sec. |3| makes it clear that the term with $ 2 in eq. ( 3.17| ) 
should be evaluated at a rapidity scale that is the earlier of the two rapidity scales y p and 
y q in the evolution of the corresponding nucleus. This prescription guarantees that the 
same is true when the scale is parametrized in terms of x instead of rapidity. 

The solution of the BK equation is reliable when the gluon density is large. The initial 
condition for the evolution is typically set at x < 0.01. For larger values of x, one expects 
the BK description to break down; we use instead a phenomenological extrapolation (used 
previously in [47, 66] ) for the unintegrated gluon distribution which has the form 

1-*^ 



cj>(x,k ± ) = 0(x o ,k ± ), (5.2) 

V 1 - / 

where xq = 0.01 and the parameter [3 = 4. This extrapolation to large x is unreliable 
and depends on physics which is not amenable to the renormalization group approach 
advocated here. However, in experiments with finite kinematic reach, it is inevitable that 



13 Another option would be to replace the momentum in eq. (5.1) by Q B (x). We have tried this and found 
our results to be insensitive to the choice of scale. 
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one is sensitive to the non-perturbative physics at large x in some kinematic range. For 
example, from the kinematic expressions in eq. ( |5,1| ), the unintegrated gluon distribution 
of gluons having p± = 0.5 GeV at RHIC energies of y/s = 200 GeV/nucleon will begin to 
be sensitive to the large x extrapolation of the distribution at y p ~ 1.4 units in rapidity. At 
the LHC energy, the range in rapidity where we avoid this sensitivity is much greater. At 

= 5.5 TeV, the same gluon does not probe the large x extrapolation of the unintegrated 
gluon distribution until y p = 4.7 

With these caveats in mind, we shall now examine the two gluon inclusive distributions 
in A+A collisions both at RHIC (y/s = 200 GeV) and at the LHC (s/s = 5.5 TeV). The 
beam rapidities for these energies are Ibeam ~ dbln ( M — ] ~ 5.36 and 8.68 for RHIC 

Y lVl nucleoli J 

and LHC respectively. We will first consider RHIC collisions and compare our results to 
recently measured long range rapidity correlations in the near-side ridge by the PHOBOS 
collaboration. The experimental quantity of interest is g5n> where the trigger particle 
consists of all particles having p± > 2.5 GeV and an acceptance in rapidity in the range 
< ^ tri s- < 1.5. The particles associated with this trigger have momenta larger than 4 
(35) MeV at a rapidity of 3 (0). In performing the At/ projection in the experiment, the 
near side yield is integrated over |A0 p(? | < 1. Hence in computing the per-trigger yield, we 
should in principle also integrate our two particle correlation C(p,q) over the PHOBOS 
acceptance. We will instead perform a more qualitative comparison here by computing 
instead our two particle correlation at representative values of the trigger and associated 
particle momenta and multiplying the result by the phase space volume corresponding to 
the PHOBOS acceptance. 




-6 -5 -4 -3-2-1012 



Figure 6: Comparison of our results for long range rapidity correlations to data from the PHOBOS 



collaboration [19]. The curves shown are obtained by adding our result (expressed by eq. (5.3) for 
long range rapidity correlations in the PHOBOS acceptance to the short range jet correlation in 
p+p collisions obtained using PYTHIA. 

For the trigger particle we take p± = 2.5 GeV at ytrig. = 0,0.75, and 1.5 units in 
rapidity. We assume the associated particle has mean p± = 350 MeV. For all cases, we 
compute the yield at A4> pq = 0. Then, in terms of our expression for the two particle 
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cumulant, the required quantity can be written as 

1 diV n ,C(p t _ L rig -,pr° c -,y trig .,yassoc. =ytrig.+A7/,A^ = 0) 



(5.3) 

where C(p, q) is the two particle cumulant given by eq. ( 3.17| ). In the above expression, the 
phase space volume corresponding to the trigger particle cancels out; we are left with an 
overall factor from the associated particle's phase space volume, Vp| soc ' , which we estimate 
to be Vp| soc ' = 7r GeV 2 . We arrive at this estimate by performing the angular integration 
over (/'assoc. times the p± integration over the acceptance. Other than Vp^ soc ', the only 
additional parameter in our expression is a s (Qg) which we take to be a s = 0.35. With 
these stated values of a s and Vp| soc ', our overall normalization is now fixed. The function 
F(A<j) pq ) comes from the collimation of the Glasma flux tubes due to radial flow as discussed 
in Ref. [34]. At A<p pq = 0, this can be expressed as 

F(A<p pq = 0) = cosh^anh" 1 0] (5.4) 

where (3 = V/c is the radial flow velocity. 

To take into account the short range correlation from fragmentation not included in 
our formalism we add to eq. (|5.3| ) the short range jet correlation resulting from PYTHIA. 
The result is compared to the PHOBOS experimental data [19] in fig. |(| One can see 
that the agreement with data is quite good. In principle the collimation from radial flow 
through eq. (|5.4 ) can be a function of rapidity. We have estimated this effect by assuming 
that the space-time and momentum space rapidity are strongly correlated. From fits to 
BRAHMS data [67-69] on the inclusive hadron spectrum, we estimate the r\ dependence 
of the flow velocity to be P(rj) = 0.72 — 0.04|^|. When including this rapidity dependent 
flow through eq. ( |5.4j ), the effect is so small that it would not result in a visible change to 
the curves plotted in figure ||[ 

At RHIC energies, the range in rapidity where the results are sensitive to small x 
physics exclusively is quite limited. At the LHC, this range is much larger and the effects 
of QCD evolution on long range rapidity correlations is more transparent. In fig. |7], we show 
results for the two particle cumulant C(p, q) as a function of the rapidity difference between 
the two gluons. In the figure on the left, the correlation is plotted for p± = q± = 2 GeV; 
the right figure corresponds to the asymmetric case of p± = 10 GeV and q± = 2 GeV. 
For both scenarios, we show the evolution in rapidity of the two particle correlation at 
different trigger rapidities y p . The solid part of each curve corresponds to the kinematic 
range where only x < 0.01 values in each of the nuclear wavefunctions are being probed. In 
contrast, the dashed part of each curve denotes the kinematic range which is sensitive to 
x > 0.01 for at least one of the nuclei; in this regime, the results are more sensitive to the 
form chosen for the large x extrapolation than to the high energy QCD evolution equations 
at small x. Because of the large kinematic reach of the LHC, we observe in fig. [7] that we 
have a region contributing to the double inclusive rapidity spectrum, of nearly 7 units in 
the rapidity difference of the two gluons, which is sensitive only to the small x evolution in 
the nuclear wavefunctions. The shape and magnitude of these correlations will therefore 
give us unique insight into the evolution of multi-parton correlations in high energy QCD. 
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Figure 7: The predicted two particle correlation spectrum as a function of the rapidity difference 
between the two gluons. The figure on the left corresponds to the case where the transverse momenta 
of the two gluons are equal and are p± , q± = 2 GeV; the figure on the right depicts the case where 
px = 10 GeV and q± = 2 GeV. The different plots reflect different trigger rapidities. Solid parts of 
each curve correspond to x < 0.01 in both nuclei; the dashed parts are sensitive to x > 0.01 in at 
least one of the nuclei. We have rescaled some curves, by the given factors, for clarity. 



6. Summary 

In Ref. [1], a general formula (eq. ( j2.6|) ) was derived for double inclusive gluon production 
in the Glasma at arbitrary rapidity separations. In this paper, we showed that this formula 
reduces to a compact expression, eq. ( 3.1 7| ) , in terms of the unintegrated gluon distributions 



in the two nuclei. This simplification holds when p±, q± > Q s and when the mean field 
Balitsky-Kovchegov (BK) framework — valid in the large iV c limit — is used to describe 
the high energy evolution of the nuclei. The unintegrated gluon distributions at small x 
are simply related to the dipole foward amplitude, which in turn satisfies the BK equation. 
Solving the running coupling form of the BK equation, with initial conditions determined 
from fits to fixed target e+A data, we computed the double inclusive spectrum at RHIC and 
LHC energies. In the case of the former, we obtained a good agreement with the PHOBOS 
data, albeit the kinematic region where small x partons are probed in both nuclei is rather 
small. In the latter case, we showed that there is a wide kinematic window for rapidity 
correlations at the LHC. Our results therefore open a new window into the study of the 
high energy evolution of multiparton correlations in nuclear wavefunctions. 
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A. Evaluation of eq. ( |3.9|) 



In this appendix, we work out some of the details of the derivation of the double inclusive 
spectrum. In particular, eq. (|3.9|) is expressed as the product of eight color charge densities. 
For a non-local Gaussian distribution of these sources, one has nine possible pairings of 
color source densities. These are evaluated explicitly below. 

^) = (27r) 8 <W^<W*<5(k 1± - k 2± )fi 2 Ai (y p , |k 1± |)<5(k 3± - k 4± )/4(y ? , |k 3 ±|) 
x 5(k 2± - k x _L)/i| 2 (y p , |p_L - ku_|) 

x 5(k 3± - k 4± )fi 2 A2 (y q , |q ± - k 3 _L|) , (A.l) 



= (2TT) s 6 b f6 dh 5 ce 6^6(k 1± - k 2± )v 2 Al (y p , |k 1± |)«5(k 3i . - k A± )^ 2 Al (y q , \k 3± \) 
x <Hq± +P± -k 3 i - ki_L)/x^ 2 (y 9 , |p ± -kul) 

x <Hq± +P± -k 4 ± - k 2 ±)fj 2 A2 (y q , |p_l — k 2J _ | ) , (A.2) 

= ( 2 7r) 8 <WVW(k 1± - k 2± )^ 2 Al (y p , \k 1± \)S(k 3± - k A± )^ 2 Al (y q , |k 3± |) 
x S{p± - q± - k 2± + k 3i _)^ 2 M (y q , |p ± - k 2 _L I) 

x 5(q ± - p_L - k 4± + ki_L)/4 2 (y 9 , |p ± - ku_|) , (A.3) 



= (2<K) 8 5f h 5 bd 5 c V5 ie 5(k 2± + k 4± )^ 2 Ai (y p , \k 2± \)5(k 1± + k 3± )^ 2 Ai (y p , \k 1± \) 
x 5(k 1± - k 2± )fi 2 A2 (y p , |p ± - ki_L |) 

x 5(k 3± - k A1 _)ii 2 M (y q , |q± - k 3± |) , (A.4) 

= (27r) 8 <WW e <Kk 1± - k 4± )fi 2 Al (y p , \k 1± \)5(k 2± - k 3± )» 2 Al (y p , |k 2± |) 
x 5(k 1± - k 21 _)n\ 2 (y p , |p ± - ki_L |) 

x <5(k 3 _L - k 4± )fi 2 A2 (y q , | q_L — k 3 _L | ) , (A.5) 

= (27r) s 5^5 bd S^ c 5(k 2± + k 4± )» 2 Ai (y p , \k 2± \)5(k 1± + k 3± )^ 2 Ai (y p , \k 1± \) 
x 5(p ± - q± - k 2 _L + k 3A _)n 2 M (y q , |p ± - k 2 _L I) 

x 5(q ± - p± - k 4± + ki_L)/4 2 (y 9 , |p± - ki_L |) , (A.6) 

= (27r) 8 <W^ c W(k 2± - k 3± )/4(y p , |k 2± |)ff(k 1± - k 4± )^ Al (y P , \k 1± \) 
x 5(q_L + p_L - k 3± - k 1± )^ A2 (y q , \p ± - k 1± \) 

x 5(q ± + p ± - k 4± - k 2 ±)/JL 2 A2 (y q , \p± - k 2± \) , (A.7) 

= (2n) s 5f h 5 bd 6 ce 5^5(k 2± + k 4± )^(y p , \k 2± \)5(k 1± + k 3± )fi 2 Ai (y p , \k 1± \) 
x 5(q± + p_L - k 3± - k 1± )/jL A2 (y q , \p± - k 1± \) 

x 5(q± + p ± - k 4± - k 2± )fjL 2 A2 (y q , \p±-k 2± \) , (A.8) 
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and 



^ = (27r) 8 5^5 bh 5^ c 5(k 2± - k 3± )/4(y p , \k 2± \)5(k 1± - k 4A _)(j, Al (y p , |k 1± |) 
x $(P± - q± - k 2± + k 3 ±)/xi 2 (y g , | pi. - k 2 i|) 

x <5(<U - P± - k 4 i + ki ± )^ s (y 9 , |pi - k 1± |) . (A.9) 

The classification of these contributions was examined previously in [34] in the frame- 
work of the MV model. The analysis is identical here. The expression is trivial as 
it cancels the square of the single particle distribution. Let us look at the <5-functions 
in jf( 4 )'( 8 ) . These yield a local [£(Pi =L qi)] 2 -contribution that we shall neglect here as 
in Ref. [34]. Similarly, expressions are sub-dominant 14 . The leading terms are 

therefore ^ r ( 1 )>( 2 )'( 3 )'( 6 ) . If we plug these back into eq. ( |3.11| ), we obtain the following four 
contributions to the two gluon spectrum 

x A t i 1 (yp>l k i_ L |)//^ s (?/g,|pi. -kii|) n 2 Ai (y q ,\p± + q ± - k 1± \) , (A.10) 
where {k^ } = {k 1± , k 1± , p ± + q ± - k 1± , p ± + q x - k^}, 

^(P,q) = ^e(^-l)5i J ^ g (pjq;{k (3) }) 

x VaAVp, |ki±|)/u^ 2 (y p , | p ± - kii|) /J,%(y q , |q± + ku_|) , (A.ll) 

(3) 

where {k\/} = {k 1± , k 1± , -kii, -kii}, 

C ,.) M .«!t^| dV6(M ; {k S' } ) 

x lA^Vpi |kii|)/4 2 (y g , |Pi -kii|) MAi^gJqi ~ P± + k ii|) , (A.12) 

(2) 

where {k^} = {k 1± , k 1± , q ± - p ± + k 1± , q ± - p ± + k 1± }, and finally 

cW(p>q) = 9 12 NHN^-1)S ± J d2ki±g(p)q;{k (6) }) 

x ^Ai(2/p 5 |kii|)^ 2 (y p , | p ± - kii|) n 2 M {y q , |q± - kij_|) , (A.13) 

where {k^} = {kii, kii, kii, kii}. Using eq. (fT|), we can express 
as eq. (gT?]). 

14 We note that there is an order one contribution coming from .T^ 5 '^ 7 ) when the relative angle between 
p±,q± is A<j> pq < Sa.. In the limit where ^J- C 1 these contributions will be washed out by re-scattering 
in the same manner as the 5-function contributions coming from .T^ 4 )*' 8 ^ and thereby not alter our result. 
We thank Kirill Tuchin for pointing out this subtlety to us. 
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B. Initial conditions for BK evolution 



The initial conditions for BK evolution of protons and nuclei are obtained by comparing 
results for the dipole cross-section to deep inelastic scattering data. The inclusive structure 
function F2 is given by 



Ff(x,Q 2 ) 



Q 2 



4vr 2 



a,- 



(B.l) 



where a A ' is the virtual photon- nucleus cross section for transverse and longitudinal po- 
larizations of the virtual photon. These in turn are given by 

a^ L (x,Q 2 ) = j dz f d 2 bd 2 r \^ L (z,Q 2 ,r)\ 2 AT A (b,r, x) , (B.2) 

where M A is the dipole-nucleus scattering amplitude. We assume here that the b depen- 
dence can be factorized as 



M A (h,v,x)=2T A (h)M A (v,x) 
The virtual photon-nucleus cross section can then be expressed as 



(B.3) 



a T A ' L { Xl Q 2 ) = a A j\z J d 2 r|* T , L (z,QV)| 2 A/- A (r,x) (B.4) 



where 



a A = 2j d 2 bf A (b) 



(B.5) 



Initial condition for protons 



The initial condition for protons was determined from a global fit of Fi data in the work 
of [43]. Two different models for the initial condition were used in that work. The first is 
the GBW model 



N(r, Y = 0) = 1 - exp 
and the other is the MV model 



2 ^2-1 



(B.6) 



N(r, Y = 0) = 1 - exp 



In 



1 



J 2K QCD 



+ e 



(B.7) 



The fit parameters obtained in [43] are summarized in the table [l]. 
Initial condition for nuclei 

We shall now consider the initial conditions for nuclei using the same model as the initial 
conditions for protons. We do not attempt to perform a global fit since the data for DIS 
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i.e. 


a p (fm 2 ) 


Q 2 so, v (GeV 2 ) 


C 2 


7 


GBW 


3.159 


0.24 


5.3 


NA 


MV 


3.277 


0.15 


6.5 


1.13 



Table 1: Parameters for the initial condition of the proton dipole cross section obtained in [43]. 




Figure 8: DIS fixed target e+A data on the ratio of structure functions as a function of A for 
fixed x = 0.0125. The plots correspond to (right) MV model initial condition having 7 = 1 and 
(left) GBW initial condition. 



off nuclei are not nearly of the same quality. We use a model where the initial saturation 
scale scales linearly with A 1 / 3 , 

Q 2 S0 =cAV 3 Q 2 S0tP , (B.8) 

where c is a constant to be determined from the data. 

In order to constrain the initial condition, we begin by looking at the New Muon 
Collaboration's (NMC) data [64] for F^/Ff as a function of A at x = 0.0125 which is 
close to our xq = 0.01. In this case, there is no BK evolution, and we have a direct 
comparison of the nuclei's initial condition with the data. In computing F^/F^ we will 

. 2/3 n-i 

need a model for how the cross section scales with A. We take a A = f-f|) a c . Figure || 
shows the NMC data as a function of A for the GBW initial condition for four different 
values of c (left) and the MV model initial condition having anomalous dimension 7 = 1 
(right). It is clear that in order to be consistent with the data we must take c ~ 0.25. The 
nuclear saturation scale given by eq. ( |B.8| ) is too small to be consistent with measurements 
by other groups. 

In fig. ||, we show the NMC data on the ratio of structure functions as a function 
of A now using the MV initial condition with anomalous dimension 7 = 1.13. We find 
that for c = 0.5 this fits the data rather well. Based on the above results, for nuclei we 
will use the MV model and take 7 = 1.13 and Q 2 = 0.5^ 1/3 Q 2 o p . Therefore we have 
Q 2 sq = 0.17, 0.26, 0.37 and 0.44 (GeV) 2 for C, Ca, Sn and Au respectively. Note that these 
values of the saturation scale are for quarks in the fundamental representation. For gluons 
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X 




Figure 9: Left: The saturation scale Q s as a function of the square root of the rapidity y 1 / 2 for 
protons, calcium and gold nuclei. Note that the slopes of the three curves approach the same value 
at large Y. Right: A = dlnQ 2 /dY as a function of Y approaches a universal value at large Y. 

in the adjoint representation, the corresponding saturation scale is 

(Qs 2 ) g = ^(Qs 2 ) q = 2.25(Q^) q (B.9) 

F 

For gold nuclei this yields (Q s 2 )g ~ 1 GeV 2 at x = 0.01, in fairly close agreement to 
the value of 1.3 GeV 2 obtained in [18,70]. Finally, we plot in fig. || (left) the saturation 
scale in the running coupling function of Y 1 / 2 for the proton, calcium and gold 

nuclei. The behavior at small Y 1 / 2 is sensitive to the initial conditions of each of these 
nuclei; however, at large Y 1//2 (small x) the curves of the three nuclear approach the same 
slope, as one expects asymptotically for the behavior of Q s when running coupling effects 
are accounted for (see also [71]). The same trend can be observed by plotting (see fig. ^ 
right) A = dlnQ 2 /dY, the parameter that sets the rate at which the dipole amplitudes 
evolve with rapidity. These results confirm the universal behavior at large Y predicted in 
Ref. [72]. 
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